/**
 * @file Interpolation.h
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2020-10-18
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _INTERPOLATIONOPERATER_H_
#define _INTERPOLATIONOPERATER_H_

#include <iostream>
#include <cmath>
#include <vector>
namespace MyMultigrid
{
    class InterpolationOperator
    {
    public:
        std::vector<double> input;
    public:
        InterpolationOperator(/* args */){};
        //InterpolationOperator::InterpolationOperator(std::vector<double> v){input = v;};
        ~InterpolationOperator(){};
        virtual void do_Interpolation(int l) = 0;
        virtual void set(std::vector<double> v) = 0;
    };
    
    class LinearInterpolaionOperator:public InterpolationOperator
    {
    public:
        LinearInterpolaionOperator(/* args */){};
        LinearInterpolaionOperator(std::vector<double> v){LinearInterpolaionOperator::input = v;};
        ~LinearInterpolaionOperator(){};
        void do_Interpolation(int l)
        {
            int low_h = pow(2,l + 1) + l - 1;
            int high_h = pow(2,l + 2) + l;
            int low_2h = pow(2,l) + l - 2;
            int high_2h = low_h - 1;
            input[high_h] = input[high_2h];
            input[low_h] = input[low_2h];
            for(int i = low_h + 1;i < high_h;i++)
            {
                int j = i - low_h;
                if(j % 2 == 0)
                    input[i] = input[low_2h + j / 2];
                else
                    input[i] = 0.5 *(input[low_2h + j / 2] + input[low_2h + j / 2 + 1]);
                
            }
        };
        void set(std::vector<double> v){LinearInterpolaionOperator::input = v;};
    };

    class QuadraticInterpolaionOperator:public InterpolationOperator
    {
    public:  
        QuadraticInterpolaionOperator(/* args */){};
        QuadraticInterpolaionOperator(std::vector<double> v){QuadraticInterpolaionOperator::input = v;};
        ~QuadraticInterpolaionOperator(){};
        void do_Interpolation(int l)
        {
           
            int low_h = pow(2,l + 1) + l - 1;
            int high_h = pow(2,l + 2) + l - 1;
            int low_2h = pow(2,l) + l - 2;
            int high_2h = low_h - 1; 
            std::vector<bool> visited(input.size(),0);
            input[high_h] = input[high_2h];
            input[low_h] = input[low_2h];
            visited[low_h] = visited[high_h] = 1;
            double h_2 = 1 / double(pow(2,l));
            double h_1 = h_2 / 2.0; 
            int i = low_h + 1;
            while (i < high_h)
            {
                int j = i - low_h;
                if(visited[i])
                    i++;
                else
                {
                    if(j % 2 == 0)
                    {
                        input[i] = input[low_2h + j / 2];
                        visited[i] = 1;
                        i ++;
                    }
                     else
                    {
                        if(i != high_h - 1)
                        {
                            int k = i + 2;
                            int x0 = low_2h + j / 2;
                            int x1 = low_2h + j / 2 + 1;
                            int x2 = low_2h + j / 2 + 2;
                            double x0_p = 0.0 + h_2 *(x0 - low_2h);
                            double x1_p = 0.0 + h_2 *(x1 - low_2h);
                            double x2_p = 0.0 + h_2 *(x2 - low_2h);
                            double y0 = input[x0];
                            double y1 = input[x1];
                            double y2 = input[x2];
                            double i_p = 0.0 + j * h_1;
                            double k_p = 0.0 + (j + 2) * h_1;
                            input[i] = y0 *(i_p - x1_p) * (i_p - x2_p) / (double(x0_p - x1_p) * (x0_p - x2_p));
                            input[i] +=y1 *(i_p - x0_p) * (i_p - x2_p) / (double(x1_p - x0_p) * (x1_p - x2_p));
                            input[i] +=y2 *(i_p - x0_p) * (i_p - x1_p) / (double(x2_p - x0_p) * (x2_p - x1_p));
                            input[k] = y0 *(k_p - x1_p) * (k_p - x2_p) / (double(x0_p - x1_p) * (x0_p - x2_p));
                            input[k] +=y1 *(k_p - x0_p) * (k_p - x2_p) / (double(x1_p - x0_p) * (x1_p - x2_p));
                            input[k] +=y2 *(k_p - x0_p) * (k_p - x1_p) / (double(x2_p - x0_p) * (x2_p - x1_p));
                            visited[i] = visited[k] = 1;
                            i++;
                        }
                    }
                }
            }
        }
         void set(std::vector<double> v){QuadraticInterpolaionOperator::input = v;};
    };
    
    
} // namespace MyMultigrid

#else

#endif
